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Abstract 

An epidemic multi-strain model with temporary cross-immunity shows 
chaos, even in a previously unexpected parameter region. Especially dengue 
fever models with strong enhanced infectivity on secondary infection have 
previously shown deterministic chaos motivated by experimental findings 
of antibody-dependent-enhancement (ADE). Including temporary cross- 
immunity in such models, which is common knowledge among field re- 
searchers in dengue, we find a deterministically chaotic attractor in the 
more realistic parameter region of reduced infectivity on secondary infec- 
tion ("inverse ADE" parameter region). This is realistic for dengue fever 
since on second infection people are more likely to be hospitalized, hence 
do not contribute to the force of infection as much as people with first 
infection. 

Our finding has wider implications beyond dengue in any multi-strain 
epidemiological systems with altered infectivity uppon secondary infection, 
since we can relax the condition of rather high infectivity on secondary in- 
fection previously required for deterministic chaos. For dengue the finding 
of wide ranges of chaotic attractors open new ways to analysis of existing 
data sets. 

1 Introduction 

We observe deterministically chaotic attractors [1, 3, 2] like the one in Fig. 1 
for a multi-strain model with less infectivity for secondary infection as for the 
first just by adding temporary cross-immunity to previously existing dengue 
models. 
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Either researchers previously focussed in their models on higher infectivity 
for secondary infection as for the first due to the hypothesized effect of antibody- 
dependent-enhancement (ADE), wich is confirmed in tissue experiments [4, 5], 
to increase viral load on a secondary infection with a different strain than 
obtained in the first infection [6, 7]. Or they focussed on temporary cross- 
immunity versus ADE but again limiting the effect of ADE to increase the 
contribution of secondary cases to the force of infection [8] . 
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Figure 1: Deterministically chaotic attractor obtained from a two-strain model 
with temporary cross-immunity in the "inverse ADE" parameter region of con- 
tribution of secondary infection to the overall force of infection. A further de- 
scription of this graphics will be given in the next sections. 

Our model is a basic two-strain SIR-type model for the host population only 
slightly refined as opposed to previously suggested models for dengue fever. It is 
capturing the effective dynamics of the human host population for the dengue 
virus, keeping effects of the mosquitoe dynamics (dengue is transmitted by 
mosquitoes acting as transmission vectors for the virus) or seasonality only in 
account by the effective parameters in the SIR-type model, but not modelling 
these mechanisms explicitly. Instead we focus on the multi-strain aspect and 
its effects on the host population. 

This basic model structure allows to generalize our findings to other multi- 
strain epidemiological systems, but is motivated by modelling dengue fever 
epidemiology with its peculiar phenomenology. 

Dengue is a mosquito-borne infection which in recent years has become a 
major international public health concern. Two fifths of the world's population 
is at risk from dengue. The World Health Organization (WHO) estimates there 
may be 50 million cases of dengue infection worldwide every year. The disease is 
now endemic in over 100 countries throughout the Americas, South-East Asia, 
the western Pacific islands, Africa and the Eastern Mediterranean [9]. 

Dengue fever is transmitted by the female domestic mosquito Aedes aegypti, 
although At. albopictus and Ae. polynesiensis can also act as vector [10]. Virus 
transmission in its simplest form involves the ingestion of viremic blood by 
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mosquitoes and passage to a second susceptible human host. The mosquito 
becomes infected when taking a blood meal from a viraemic person. After 
an extrinsic incubation period, the mosquito becomes infective and remains so 
during its entire life span [12]. 

As the blood meal stimulates ovoposition, which undergoes at least one, 
often more, reproductive cycles there is an opportunity of vertical transmission 
to the eggs, passing the virus to the next generation of mosquitoes [13, 11]. In 
a dry state, eggs survive for very long periods and the virus can persist until 
the subsequent rainy season [14]. 

Dengue is caused by four antigenically distinct single-strand positive-polarity 
RNA viruses, designated dengue types 1, 2, 3, and 4, and belonging to the fam- 
ily Flaviviridae [9]. Infection by one serotype confers life- long immunity to 
only that serotype and temporary cross-immunity to other serotypes exists. It 
lasts from three to nine months, when the antibody levels created during the 
response to that infection would be enough to protect against infection by a 
different but related serotype [17, 8]. 

Two forms of the disease exist: dengue fever (DF) or classic dengue, of- 
ten benign, and dengue hemorrhagic fever (DHF), which may evolve towards 
a severe form known as dengue shock syndrome (DSS) [15]. Without proper 
treatment DHF case fatality rates can exceed 20% [9]. DF is characterized by 
headache, retro-orbital pain, myalgia, arthralgia, rash, leukopenia, and mild 
thrombocytopenia. The symptoms resolve after 2-7 days. DHF is a poten- 
tially deadly complication that is characterized by high fever and haemorrhagic 
phenomenae. DHF develops rapidly, usually over a period of hours, and re- 
solves within 1-2 days in patients who receive appropriate fluid resuscitation. 
Otherwise, it can quickly progress to shock [16]. 

There are indeed pre-existing antibodies to previous dengue virus that can- 
not neutralize but rather enhance infection in vitro, a process described as 
antibody-dependent enhancement. Epidemiological studies support the associ- 
ation of DHF and DSS with secondary dengue infection. Halstead [9, 4, 5] found 
that DHF and DSS were 15-80 times more likely in secondary then in primary 
infections and were positively associated with pre-existing dengue- virus-specific 
antibodies. However, there is no animal model of DHF and DSS, and the the 
causal relationship between ADE and severe disease remains unverified [16]. 

There is no specific treatment for dengue. A vaccine against dengue is not 
yet available, since it would have to simulate a protective immune response to 
all four serotypes [18]. 

Mathematical models describing the transmission of dengue viruses ap- 
peared in the literature as early as 1970 [19]. More recently modelling attention 
has focussed on including ADE as due to higher viral load of hosts on secondary 
infection than on the first, hence a higher contribution to the force of infection 
of each strain, reporting deterministically chaotic attractors [6] and chaos de- 
synchronization [7] to explain the co-existence of the known four dengue viral 
strains. Temporary cross-immunity against all strains after a first infection has 
been included in mathematical models as well [8], but to our knowledge, no 
systematic investigation of the attractor structures of simple two-strain models 
with dengue-realistic temporary cross-immunity and decreased contribution of 
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secondary infection to the force of infection, due to severity of infection with a 
second strain and eventual hospitalization, has been performed so far. 

Our finding of chaotic attractors in this parameter region thus indicate that 
deterministic chaos is wider present than previously expected in multi-strain 
models. For such scenarios new tools of non-linear data analysis like Takens' 
embedding are available [20, 21], and allow to obtain topological information 
(fixed points, periodic orbits and the nature of chaotic attractors) about the 
whole multi-strain epidemiological system from time series of overall infecteds 
only, not needing any single strain data sets. 

In the next section we present the two-strain model and give the deteministic 
mean field version of it. We then analyze the deterministic ODE-model in 
various parameter regions in its attractor structure. 



2 The basic epidemic model 

Multi-strain dynamics are modelled with SIR-type models, dividing the host 
population into susceptible, infected and recovered individuals. 

In the simple SIR epidemics without strain structure of the pathogens we 
have the following reaction scheme for the possible transitions from one to 
another state 

s + i -!U i + i 

I R 

R S 
S, I, R — > S 

for a host population of N individuals, with contact and infection rate (3, re- 
covery rate 7 and temporary immunity rate a. Demography is denoted in the 
last reaction as exits from all classes S, I and R with rate n to the new born 
susceptibles. 

To include population noise stochastic models are investigated. For state 
vectors n, here for the SIR- model n = (S, I, R), the master equation [22] reads 

with transition rates Wn t n given by the following expressions 

W(S+1,I,R-1),(S,I,R) = ol-R 

W (S -1,I+1,R),(S,I,R) = P ■ ^ S ( 2 ) 
W (S,I-1,R+1),(S,I,R) = I' 1 
from which the rates Wn,n follow immediately as 

W(s,i,r),(s-i,i,r+i) = a-(R + l) 

w {s ,i,rus+i,i-i,r) = (3 ■ (S + 1) (3) 

W(S,I,R),(S,I+l,R-l) = 1 ■ (I + 1) • 
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In addition we have the transitions for the demography with rate fi, e.g. for 
the transition from recovered to susceptibles wts+i,i,R-i),(S,i,R) = ft ■ R- This 
formulation defines the stochastic process completely to capture demographic 
noise and will be the basis for the multi-strain model to be described in the 
following. 

The above mentioned deterministic ODE model describes in mean field ap- 
proximation the dynamics of the mean values, e.g. 

N 

(I):=J2IP(SJ,R) ■ (4) 

7=0 

The dynamics for the mean value is then given by inserting the master equation 
j t {I) = £ ■ (I-S) - ( 7 + „) ■ (/) » | ■ </) ■ (S) - ( 7 + ,x) ■ (/) (5) 

with the mean field approximation (I ■ S) ~ (/} • (S). Carried out for all mean 
values (S), (I) and (R) gives a closed ODE system for the SIR dynamics, which 
is a deterministic system in the sense that initial values determin the time course 
of the system for all times. The mean field ODE system reads, now omitting 
the brackets for mean values, 

^ = a R- — -I-S + fi(N-S) 
dt N 

dR T 

—— = 7-T — aR — aR 
at 

and has only fixed points as possible stationary solutions, i.e. attractors. Tran- 
sients under certain parameter values oscillate into the fixed point, hence can 
be already more complex than the final attractor. Stochastic versions of such 
models with only fixed points possible as attractors but oscillating transients 
are reported to also show stabilization of the oscillations due to population noise 
[23, 24]. So the deterministic mean field ODE models with their attractors give 
rather the minimal complexity caused by the model structure. 



2.1 Basic two-strain model 

For two strains, 1 and 2, we have the following SIR-type model, now labelling 
the SIR classes for the hosts that have seen the individual strains. Susceptibles 
to both strains (S) get infected with strain 1 or strain 2 (I2), with force 
of infection [3\ and fa respectively. They recover from infection with strain 
1 (becoming R\) or from strain 2 (becoming R2), with recovery rate 7. In 
this recovered class, people have full immunity against the strain that they 
were exposed to and infected, and also, temporary immunity against the other 
strain (called period of temporary cross-immunity). After this, with rate a, 
they enter again in the susceptible classes (Si respectively S2), where the index 
represents the first infective strain. Now, Si can be reinfected with strain 2 
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(I12) and S2 can be reinfected with strain 1 (I12), with infection rates fafo 
and (f>iPi- The parameter <p in our model acts as "inverse ADE parameter", 
decreasing the infectivity of secondary infection, where people are more likely to 
be hospitalized because of the severity of the disease (DHF), do not contributing 
to the force of infection as much as people with first infection do. Finally, I12 
and I\2 go to the recovered class (R), immune against all strains. We include 
demography of the host population by denoting the birth and death rate by /i, 
assuming constant population size N. 

In the following reaction scheme we describe the transitions for first infection 
with strain 1 and secondary infection with strain 2 

h + h 
h + hi 
Ri 
Si 

In + h 

I12 + 112 
R 

and for the reverse process, where the first infection is caused by strain 2 and 
the secondary infection is caused by strain 1. 

S + h 

S + I12 
h 
R2 
S2 + I1 
S2 + hi 
hi 

The demographic transitions are 

S, h, h, Ri, R2, Si, S2, 112, hi, R S 

defining the system of two strains completely. 

The stochastic version is now in complete analogy to the previous SIR model, 
and we can calculate the mean field equations (or simply read it off from the 
reaction scheme directly). 



S + h 

S + hi 
h 
Ri 
Si+h 
Si + 112 

I12 



Pi 

7 ; 

a 
7 



02/32 
7 ; 
a 

Pi , 
<hfh 
7 



h + h 
h + 112 
R2 
S 2 

hi + h 

hi + hi 
R . 
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The mean field ODE system for the two strain epidemiological system is 



§ = -§S(/i + ^i/2i)-§5(7 2 + ^/ 12 )+ M (iV-5) 
^ = |s(/ 1 + ^/ 21 )- ( 7 + n)h 
^ = ^5(J 2 + ^/ 12 )-( 7 + A*)/2 



= 7/1 - (a + |Lt)i?i 

= 7 / 2 - (a + /x)i? 2 (7) 
Si(h + folu) + acRi - fiSi 



dRi 
~dT 

dR 2 
dt 

dSi p 2 

~dT ~ ~iv 

^ = -^S 2 (h + 4>ihi) + aR 2 - nS 2 

^ = ^5 2 (/i + ^1/21) - (7 + ^)^21 

= l{h2+hi)-^R 

We consider <f>\ = 2 = <fi, Pi = P2 = P for the moment, hence no epidemiolog- 
ical assymetry between the strains. 

To take biological information from experiences in dengue into account we 
fixe the transition rates of the model as far as is known, and only will vary 
the unknown parameters <fi and eventually a. Future work will have to also 
adjust the other parameters better to describe actual data of dengue cases. 
The parameter values are, if not otherwise explicitly stated, /j, = l/65years, 
7 = 52y~ 1 , (3 = 2 ■ 7 and a = 2y~ 1 . For the chaotic attractors we take the 
exemplaric values (/) = 0.7 and (ft = 2.7. We will also vary <j) continuously to 
obtain bifurcation diagrams [3, 2]. 



3 Analysis of the multi-strain model with temporary 
cross-immunity 

We will first look at time series simulations of the present model, Eq. system (7), 
and from there we performed a detailed analysis of chaos, investigating state 
space plots for various values for 4>, observing a rich structure of attractors 
from fixed points to bifurcating limit cycles and chaotic attractors. Finally, we 
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will show whole bifurcation diagrams for a systematic analysis for the entire 
parameter region of <f> for various values of the temporary cross-immunity a. 



3.1 Time series for </> > 1 

In order to classify the dynamic pattern of the model for various parameters 
we discard long transients which would carry some information of the initial 
conditions. So we neglect the first 2000 years in the following simulations. We 
first simulate time series for <j> > 1, see Fig. 2 and Fig. 3. 
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Figure 2: For <j) = 2.7, time series simulations of a) the suceptibles S, b) the 
at the moment infected with strain 1 (I\ + In), c) the at the moment infected 
with strain 2 (I2 + In)- Temporary cross-immunity is assumed to be a = 2. 
The absolute numbers on the y-achses indicate percentage (N = 100,). 
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Figure 3: Using the same parameter values above, In a) we plot the time series 
of the total number of infected I := I\ + I2 + I12 + I21, o,nd in b) time series of 
the logarithm of the overall infected (ln(I) ), for the same time interval. Very 
deep troughs are observed for these parameter values. 



However, also the transients reflect the dynamic behaviour of the system 
under the present parameter values. Spiraling into a fixed point often indi- 
cates a nearby periodic orbit. Irregular transients indicate chaotic behaviour 
in neighbouring parameter regions etc. Hence, even if dengue fever would be 
evolutionarily younger than the mentioned 2000 years of transient, the observed 
pattern would still give information about the dynamics of the system. Also the 
stochastic version of this deterministic mean field model would be rather more 
complex than the attractor classification and not simpler [23, 24], of course be- 
sides extinction due to population fluctuations for eventually very low numbers 
of infected. 
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Rather low troughs are observed for the total number of infected leading 
to unrealistically low numbers of infected on average in low epidemic times. 
The present deterministic model can be interpreted as a mean field model for 
stochastic models capturing the effect of population noise e.g. in terms of 
master equations as mentioned before [22]. The ODE solutions are then the 
mean values of such stochastic models. Population fluctuations would in the 
present case drive always the system to extinction. 



3.2 Time series for </> < 1 

Investigating time series for (j) < 1, as would be realistic for dengue fever due to 
more severe disease upon reinfection, hence larger chance of being hospitalized, 
also indicates more complicated dynamic behaviour than just simple fixed point 
or limit cycles. Discarding the first 2000 years, see Fig. 4 and Fig. 5. 



a) 




b) 
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Figure 4: For a = 2 and <f> = 0.7, time series simulations of a) the suceptibles 
S, b) the at the moment infected with strain 1 (I\ + hi), and c) the at the 
moment infected with strain 2 (I2 + In)- 
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Figure 5: Using the same parameter values above, in a) we plot the time series 
of the total number of infected I := I\ + 12 + I12 + I21, an d in b) time series of 
the logarithm of the overall infected (ln(I) ), for the same time interval. 



For (j) < 1 the number of infected stays quite away from zero, avoiding the 
chance of extintion in stochastic systems with reasonable system size. In Fig. 
3 b) the logarism of total infected goes as low as —70, when <f> > 1. For <fi < 1, 
see Fig. 5 b) the logarism of total infected does not pass below —7. 

This encourages to also look closer to the parameter region of < 1, the 
region of "inverse ADE", i.e. when dengue patients with severe disease because 
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of the ADE phenomenon contribute less to the force of infection due to possible 
hospitalization, and not more, as previous models suggested. 

3.3 State space plots 

In Fig. 6 we plot 500 years of dynamics to clearly observe the dynamic patterns 
present. 

41 
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a) 

Figure 6: For a = 2 and (p = 0.7, we plot in a) the time series of the susceptibles 
after discarding a long transient (here 2000 years), and in b) the time series of 
the logarithm of the overall infected, I := I\ + I2 + I12 + I21 f or the same time 
interval. In c) State space plot of the number of susceptibles versus logarithm of 
the overall infected, as show from their time series in Fig. 6 a) and b), finding 
a chaotic attractor. Parameters: a = 2, (ft = 0.7. 




We now plot the total number of infected with one strain versus the to- 
tal number of infected with the other strain. In case of synchronicity of the 
two strains we would expect the system to stay closer to the main diagonal 
of the plot, whereas for anti-synchronicity it would mainly stay near the off- 
diagonal. The observed pattern is more supporting anti-synchronicity, but near- 
synchronicity can also be observed at times (see Fig. 7). For quantification of 
these effects of chaos synchronization or chaos desynchronization applied to 
dengue models see Schwartz et al. [7]. 



-2 




12 -\ , 1 1 1 1 1 

-14 -12 -10 -8 -6 -4 -2 

Nli+W ffl 



Figure 7: State space plot of the logarithmic number of infected with strain 1 
versus the ones with strain 2. Parameters: a = 2, (f> = 0.7. 

In the further analysis we will investigate the state space plots in terms of 
the variables S and the logarithm of / := I\ + I2 + I12 + hi, since dengue 
notification data often do not distinguish between the circulating strains, hence 
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the total number of infected I, whereas the susceptible class S is N minus every 
host who ever has experienced an infection. 

In eventual data analysis the method of delay coordinates even allows to 
only work with one time series of /, and analysing I(t), I(t + r) etc. with a 
time delay r [20, 25]. 



3.4 Bifurcations of limit cycles when changing <f> 

The state space plots for various values of <fi show bifurcations from fixed point 
behaviour to limit cycles, which then bifurcate into doubledimit cycles etc., 




a) , b) 




Figure 8: Attractors for various values of (f) < 1: a) fixed point for cj) = 0.1, 
b) limit cycle for <f> = 0.2, c) bifurcating limit cycle for <\> = 0.45, d) chaotic 
attractor for <p = 0.6. 




Figure 9: Attractors for various values of 4> around 1 and larger: a) chaotic 
attractor for <j) = 0.8, b) limit cycle for <f> = 1.5, c) more complicated limit cycle 
for 4> = 1.9, d) chaotic attractor for <f> = 2.7. 
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until completely irreguar behaviour, which is the fingerprint of deterministc 
chaos. See Fig. 8 and Fig. 9. 

Looking for higher values of <fi, than just towards the first chaos window, 
shows that the chaotic attractor becomes unstable again, just leaving simple 
limit cycles as attractors for large parameter regions beyond of (p = 1. Only for 
much higher values of 4> » 1, another chaotic attractor appears, the so-called 
ADE chaotic attractor [6, 7]. 



3.5 Map of maxima of J in state space 

We also investigate maxima maps in order to classify the dynamics for various 
parameter values, see Fig. 10 and Fig. 11. 
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Figure 10: a) Limit cycle for <fi = 1.9, b) The dots indicate the local maxima 
calculated for the twisted limit cycle of a). The three dots characterize the three 
limit cycle. Parameters: a = 2, <f> = 1.9. 

We plot for the time t max , at which the total number of infected I(t) := 
h + h + In + hi has a local maximum, the number of infected at that time 
ln(I(t max )) and for the same time value the susceptibles S(t max ). In this way 
we obtain a maxima map. 
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Figure 11: In a), the determined local maxima for the chaotuc attractor found 
for a = 2 and <j> = 0.7 (See Fig. 6 c)). In b), scattered maxima plot for the 500 
years of simulation. In c) the map of maxima of the overall infected and the 
respective values for the susceptibles for a very long simulation, using the same 
conditions as for Fig 11 b). We observed that even after 200000 years, the dots 
never come back to the same point. Much more structure than the scatter plot 
appears. The fingerprint of the chaotic attractor is clearly visible now. 

Whereas in Fig. 11 b) we use 500 years of simulations, observing rather 
erratical scattering of points of local maxima, for Fig. 11 c) we use the local 
maxima of a long simulation of 200 000 years to observe the clear fingerprint of 
deterministically chaotic behaviour. 
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3.6 Bifurcation diagrams 



Ploting the local maxima of ln(I) over the varying parameter 4>, we obtain a 
bifurcation diagram, in which fixed points and simple limit cycles appear as one 
dot per parameter value, whereas double-limit cycles appear as two dots, more 
complicated limit cyles as more dots, and chaotic attractors as continuously 
distributed dots for a single 4> value (see Fig. 12). 




Figure 12: Bifurcation diagram for the maxima of the overall infected with 
changing parameter <\> and a = 2. We observed a chaotic window for <p < 
1, region of "inverse ADE", where this dynamical behaviour has never been 
described before, and another one for <fi > 1, called by ADE chaotic window, 
found in previous publications [6, 7]. 

To be sure that this unexpeted behaviour for < 1 not just appears because 
of this specific a value, we look at the robustness of the findings with varying 
the temporary cross-immunty parameter values. 

For q = 1 e.g, both chaotic windows appears, and surprisingly in the region 
of "inverse ADE" this window is even larger (see Fig. 13). 




Figure 13: Bifurcation diagram for the maxima of the overall infected with 
changing parameter <f>. Now a = 1. The bifurcation diagram appears to be quite 
robust against changes of parameters around the region under investigation. 

Only for very large values of a, where temporary cross-immunity becomes 
unimportant due to the low resident times in the classes R\ and R2, the chaos 
window for < 1 disappears, and then ADE as increasing infectivity on a 
secondary infection condition seems to be the only mechanism to observe de- 
terministic chaos (see Fig. 14). 
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Figure 14: Bifurcation diagram for the maxima of the overall infected with 
changing parameter (f>. Now for a very large value of a we get close to the models 
found in the literature, where the temporary cross-immunity is not considered, 
a) For a = 10, The chaos window for 4> < 1 disappears completely as the 
temporary cross-immunity becomes shorter or unimportant, b) For a = 20 it 
becomes clear that there is no other dynamics in the region for 4> < 1, than fixed 
points or simple limit cycles. 
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Figure 15: Bifurcation diagram for the extrema of the overall infected with 
changing parameter (p. a = 2. For the chaotic window for (f> > 1 the minimal 
values go to very low numbers of infected. 

Finally, plotting not just the local maxima but also the local minima of 
ln{I) in a bifurcation plot, as done in Fig. 15, we observed again for (p > 1 
very low troughs of infected, whereas in the chaotic region for <f> < 1 the overall 
number of infected stays always sustainably high. 

4 Summary 

Our analysis shows that also for "inverse ADE", i.e. </> < 1, deterministic chaos 
can be observed, when taking the temporary cross-immunity reported in the 
medical literature into account. 

This indicates that deterministic chaos is much more important in multi- 
strain models than previously thought. We could show this in a very basic 
model with only two strains and one reinfection possible, not needing the ADE 
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mechanism, but rather stating that upon second infection hosts spread a disease 
less likely, since it might be more harmfull. This mechanism could be present 
in other diseases than dengue fever, hence of much wider interest. 
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